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It  has  been  shown  that  the  plastic  response  of  many  materials,  including  some  metallic 
alloys,  depends  on  the  stress  state.  In  this  paper,  we  describe  a  plasticity  model  for  isotro¬ 
pic  materials,  which  is  a  function  of  the  hydrostatic  stress  as  well  as  the  second  and  third 
invariants  of  the  stress  deviator,  and  present  its  finite  element  implementation,  including 
integration  of  the  constitutive  equations  using  the  backward  Euler  method  and  formula¬ 
tion  of  the  consistent  tangent  moduli.  Special  attention  is  paid  for  the  adoption  of  the 
non-associated  flow  rule.  As  an  application,  this  model  is  calibrated  and  verified  for  a 
5083  aluminum  alloy.  Furthermore,  the  Gurson-Tvergaard-Needleman  porous  plasticity 
model,  which  is  widely  used  to  simulate  the  void  growth  process  of  ductile  fracture,  is 
extended  to  include  the  effects  of  hydrostatic  stress  and  the  third  invariant  of  stress  devi¬ 
ator  on  the  matrix  material. 

©  2010  Elsevier  Ltd.  All  rights  reserved. 


1.  Introduction 

Plasticity  describes  the  deformation  of  a  material  undergoing  non-reversible  changes  of  shape  in  response  to  applied 
forces.  Our  ancestors  in  ancient  times  already  recognized  the  plastic  behavior  of  metals  as  they  attempted  to  make  various 
tools  and  weapons.  However,  the  scientific  study  of  plasticity  may  justly  be  regarded  as  beginning  in  1864  when  Tresca  pub¬ 
lished  his  results  on  punching  and  extrusion  experiments  and  formulated  his  famous  yield  criterion  (Tresca,  1864).  This  yield 
criterion  was  then  used  by  Saint-Venant  (1870)  and  Levy  (1870)  in  their  development  of  a  theory  of  rigid-perfectly  plastic 
solid.  Another  well  known  yield  criterion  was  proposed  by  von  Mises  (1913)  on  the  basis  of  purely  mathematical  consider¬ 
ations  and  later  was  interpreted  by  Hencky  (1924)  as  plastic  yielding  occurs  when  the  elastic  shear-strain  energy  reaches  a 
critical  value.  Von  Mises  also  independently  proposed  equations  similar  to  Levy’s  for  rigid-perfectly  plastic  materials.  Other 
important  contributions  in  the  early  development  of  the  plasticity  theory  include  the  works  by  Prandtl  (1925),  Reuss  (1930), 
among  others.  Subsequently,  within  the  scope  of  elastic-plastic  materials  under  small  deformation,  the  notation  of  yield  in 
the  stress  space  formulation  was  generalized  to  cover  work-hardening  materials  and  a  unified  theory  of  plasticity  began  to 
emerge  after  World  War  II  (Hill,  1950;  Mendelson,  1968).  The  flow  theory  is  the  most  widely  known  theory  of  plasticity, 
which  consists  of  a  yield  criterion,  a  flow  rule,  a  hardening  law  and  the  loading-unloading  conditions.  The  yield  criterion 
determines  the  stress-state  when  yielding  occurs,  the  flow  rule  describes  the  increment  of  plastic  strain  after  yielding, 
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the  hardening  law  characterizes  the  evolution  of  the  flow  stress  with  increased  plastic  deformation,  and  the  loading-unload¬ 
ing  conditions  specify  if  the  stress  path  moves  outward  from  the  current  yield  surface,  inward  from  the  current  surface  or 
along  the  current  surface. 

The  most  popular  continuum  plasticity  model  is  the  so-called  J2-fl ow  theory.  This  theory  assumes  hydrostatic  stress  as 
well  as  the  third  invariant  of  the  stress  deviator  has  no  effect  on  plastic  yielding  and  the  flow  stress  and  has  been  widely 
employed  to  describe  the  plastic  response  of  metals.  However,  increasing  experimental  evidences  show  that  this  assumption 
is  invalid  for  many  materials.  Inspired  by  the  extensive  experimental  results  reported  by  Spitzig  et  al.  (1975,  1976)  on  the 
behavior  of  high-strength  metals  undergoing  uniaxial  tension  and  compression,  Brunig  (1999)  presented  an  l\-]2  yield  cri¬ 
terion,  which  is  similar  to  the  Drucker-Prager  yield  condition  in  soil  mechanics  (Drucker  and  Prager,  1952),  to  describe  the 
effect  of  the  hydrostatic  stress  on  the  plastic  flow  properties  of  metals.  Later,  the  U~]2  plasticity  model  was  applied  by  Brunig 
et  al.  (2008)  to  study  the  effect  of  stress  triaxiality  on  the  onset  and  evolution  of  damage  in  ductile  metals.  In  Brunig  et  al. 
(2000),  the  authors  added  a  J3  term  in  the  yield  function  to  study  the  deformation  and  localization  behavior  of  hydrostatic 
stress  sensitive  metals.  Kuroda  (2004)  presented  a  phenomenological  plasticity  model  accounting  for  hydrostatic  stress-sen¬ 
sitivity  and  vertex-type  of  effect.  Hu  and  Wang  (2005)  proposed  a  stress-state  dependent  yield  criterion  for  isotropic  ductile 
materials.  Cazacu  and  Barlat  (2003)  and  Soare  et  al.  (2007)  extended  Drucker’ s  J2-J3  yield  function  (Drucker,  1949)  to  include 
plastic  anisotropy  and  applied  it  to  simulate  sheet  forming.  Plunkett  et  al.  (2008)  and  Cazacu  et  al.  (2010)  extended  the  iso¬ 
tropic  yield  function  proposed  by  Cazacu  et  al.  (2006)  to  account  for  the  anisotropic  plastic  response  of  textured  metals.  Bai 
and  Wierzbicki  (2008)  discussed  a  pressure  and  Lode  dependent  metal  plasticity  model  and  its  application  in  failure  analysis. 
Gao  et  al.  (2009)  noticed  the  plastic  response  of  a  5083  aluminum  alloy  is  stress-state  dependent  and  were  forced  to  use 
different  stress-strain  curves  to  analyze  specimens  experiencing  different  stress  states.  In  a  most  recent  study,  Mirone 
and  Corallo  (2010)  found  that,  for  the  metals  they  tested,  the  hydrostatic  stress  plays  a  significant  role  in  accelerating  failure 
but  has  negligible  effect  on  the  stress-plastic  strain  relationship,  while  the  Lode  angle  has  a  considerable  effect  in  modifying 
the  stress-strain  curves  but  does  not  significantly  affect  the  failure  strains.  These  findings  agree  with  the  findings  by  Gao 
et  al.  (2009). 

A  common  approach  in  the  metal  plasticity  theory  is  the  adoption  of  the  so-called  associated  flow  rule  or  normality  con¬ 
dition,  requiring  the  yield  function  and  the  flow  potential  to  be  identical.  Although  the  normality  condition  is  a  consequence 
of  the  postulate  of  maximum  plastic  dissipation,  many  studies,  such  as  Mroz  (1963),  Nemat-Nasser  and  Shokooh  (1980),  Ne- 
mat-Nasser  (1983, 1992),  Runesson  and  Mroz  (1989),  Brunig  et  al.  (2000),  Stoughton  (2002)  and  Stoughton  and  Yoon  (2004, 
2006),  indicate  that  appropriate  constitutive  description  of  many  materials  can  be  achieved  by  using  the  less  restrictive  non- 
associated  flow  rule.  Brunig  (1999)  and  Brunig  et  al.  (2000)  demonstrated  that  pressure  sensitive  yielding  and  non-associ- 
ated  flow  rule  remarkably  influence  the  onset  of  localization  and  the  subsequent  localization  behavior.  To  model  the  ductile 
fracture  process  in  solids,  Ma  and  Kishimoto  (1998)  proposed  a  non-associated  flow  rule  to  characterize  the  yield  and  plastic 
deformation  of  void-containing  materials,  where  the  yield  function  takes  the  form  of  the  Gurson-Tvergaard-Needleman  por¬ 
ous  plasticity  model  (Gurson,  1977;  Tvergaard,  1981, 1982;  Tvergaard  and  Needleman,  1984)  while  the  flow  potential  takes  a 
slightly  different  form.  In  a  recent  publication,  Cvitanic  et  al.  (2008)  presented  detailed  finite  element  formulations  for  sheet 
metal  forming  based  on  non-associated  plasticity. 

In  this  paper,  we  describe  a  plasticity  model  for  isotropic  materials,  which  is  dependent  on  the  second  and  third  invariants 
of  the  stress  deviator  as  well  as  the  hydrostatic  stress,  and  present  its  finite  element  implementation  including  integration  of 
the  constitutive  equations  using  the  backward  Euler  method  and  the  formulation  of  the  consistent  tangent  moduli.  The  der¬ 
ivation  will  be  based  on  small-strain  formulation.  For  finite  strain  plasticity,  kinematic  transformations  are  performed  first  so 
that  the  constitutive  equations  governing  finite  deformation  are  formulated  using  strains-stresses  and  their  rates  defined  on 
an  unrotated  frame  of  reference.  Once  the  kinematic  transformations  have  eliminated  rotation  effects  on  rates  of  the  tenso- 
rial  quantities,  the  stress  updating  procedure  and  the  consistent  tangent  stiffness  formulation  remain  the  same  as  those  for 
small-strain  formulation.  Most  commercial  finite  element  programs  adopt  this  kind  of  treatment  for  finite  strain  plasticity 
thus  only  small-strain  formulation  needs  to  be  considered  in  development  of  a  user  material  subroutine.  As  an  application, 
the  proposed  plasticity  model  is  calibrated  and  verified  for  a  5083  aluminum  alloy.  Furthermore,  we  extend  the  Gurson- 
Tvergaard-Needleman  porous  plasticity  model  to  include  the  effects  of  the  stress  state  and  present  a  few  numerical  exam¬ 
ples.  The  modified  porous  plasticity  model  is  expected  to  improve  accuracy  in  predicting  ductile  fracture  process  of  certain 
materials. 


2.  Plasticity  modeling 

2.1.  Influence  of  I *  and  J3 


Let  Gij  be  the  stress  tensor  and  <7i,  g2  and  g3  be  the  principal  stress  values.  The  hydrostatic  stress  (or  mean  stress)  can  be 
expressed  as 


i 


1 

3 


(Cl  +  o2  +  o3) 


(1) 
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where  U  represents  the  first  invariant  of  the  stress  tensor  and  the  summation  convention  is  adopted  for  repeated  indices.  Let 
crjj  be  the  stress  deviator  tensor  and  o\ ,  g'2  and  o'3  be  its  principal  values,  i.e., 

a'ij  =  <7,j  -  ahSij  (2) 

where  Sy  represents  the  Kronecker  delta.  It  is  obvious  that  the  first  invariant  of  the  stress  deviator  tensor  is  zero.  The  second 
and  third  invariants  of  the  stress  deviator  tensor  are  defined  as 


h  +  OjG 3  +  <73<7'l)  =  |[(0-1  -  0-2)2  +  (<72  -  O's)2  +  (0-3  -  Cl)2] 

h  =  det(clj)  =  \o\f'jko'ki  =  a\a'2a'3 


(3) 


For  isotropic  materials,  the  plastic  behavior  is  often  described  by  the  stress  invariants  l1,J2  and  J3  and  consequently,  the 
general  forms  of  the  yield  function  (F)  and  the  flow  potential  (G)  are  expressed  as  functions  of  G,  J2  and  J3.  Eq.  (4)  describes 
the  yield  condition 


f(/iJ2J3)-c  =  o 


(4) 


where  a  is  the  hardening  parameter.  When  material  deforms  plastically,  the  inelastic  part  of  the  deformation  is  defined  by 
the  flow  rule 


,p_,dC(hj2j3) 
lJ  dav 

where  s jj.  are  the  rates  of  the  plastic  strain  components  and  A  is  a  positive  scalar  called  the  plastic  multiplier. 
A  simple  form  of  the  yield  function  is  given  as  follows 

F  =  Cl  (a,/?  +  27J| +  bj2),/6 


(5) 


(6) 


where  ai,  bA  and  Ci  are  material  constants.  The  yield  function  defined  by  Eq.  (6)  is  a  first  order  homogeneous  function  of 
stress,  from  which  an  equivalent  stress  can  be  defined  as 

Ve=F  (7) 

The  constant  cx  can  be  found  by  substituting  the  uniaxial  condition  into  Eq.  (6),  thus 

/  4  V1/6 

Ci  —  ^ai  +^29^i  +  (8) 


When  the  material  is  subjected  to  a  uniaxial  stress  cr,  the  value  of  Ci  given  by  Eq.  (8)  ensures  oe  =  o. 

The  flow  potential  takes  a  similar  form,  i.e., 

G  =  c2(a2/?  +  27j!+bz/!)1/6  (9) 

where  c2  =  (a2  +  4b2/729  +  1)-1/6  . 

If  the  flow  potential  and  the  yield  function  are  identical,  i.e.,  F=  G,  Eq.  (5)  becomes  the  so-called  associated  flow  rule.  Fur¬ 
thermore,  if  aA=b\  =  a2  =  b2  =  0,  the  plasticity  model  degenerates  to  the  formulation  of  classical  J2-flow  theory  and  ae  defined 
by  Eq.  (7)  becomes  the  von  Mises  equivalent  stress. 

The  hardening  parameter  depends  on  the  strain  history.  By  enforcing  the  equivalence  of  plastic  work,  i.e., 

oi?  =  <7,; f8?.  (10) 

the  equivalent  plastic  strain  increment  can  be  defined  as 

ep  =  (H) 

Therefore,  the  hardening  behavior  can  be  described  by  a  stress  vs.  plastic  strain  relation  g(sp),  where  sp  =  fspdt. 

Since  the  flow  potential  is  taken  to  be  a  first  order  homogeneous  function  of  stress,  Euler’s  homogeneous  function  theo¬ 
rem  results  in 


Xay—  AC  (12) 

From  (10)  and  (12),  the  plastic  multiplier  and  the  equivalent  plastic  strain  rate  can  be  related  through 

A  =  s^  =  e^  (13) 

If  the  material  follows  the  associated  flow  rule,  F=G  and  A  =  tp. 
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Because  of  the  f  term  in  the  flow  potential  G,  the  plastic  response  becomes  dilatant,  with  the  rate  of  volume  change  given 
by 

4  =  =  ^c6M/C5)  =  3  c\a2l\W  T  (14) 

2.2.  Numerical  implementation  of  the  U-J2-J3  plasticity  model 

Following  the  procedures  of  Aravas  (1987)  and  Kim  and  Gao  (2005),  we  implement  the  f  -J2-J3  plasticity  model  outlined 
in  the  previous  section  into  ABAQUS  (SIMULIA,  2008)  via  a  user  defined  subroutine. 

2.2 A.  Stress  update 

We  adopt  the  backward  Euler  method  to  integrate  the  rate  constitutive  equations.  For  the  strain  driven  integration  algo¬ 
rithm,  the  stresses  and  state  variables  are  known  at  the  start  of  each  increment  and  their  values  need  to  be  updated  at  the 
end  of  the  increment  corresponding  to  the  total  strain  increment.  The  elasticity  equations  give 

K)t+At  =  GjiMi)t+M  =  +  A skl  -  Ae£,}  =  al  -  C'jklAt?kl  (15) 

where 


^co«{(e«)t  +  A8w}  (16) 

is  the  elastic  predictor,  t  represents  the  time  at  the  start  of  the  increment,  t  +  At  represents  the  time  at  the  end  of  the  incre¬ 
ment,  and  the  superscripts  e  and  p  denote  elastic  and  plastic  components,  respectively.  The  total  strain  increment  A ew  is 
known,  and  if  the  linear  elastic  behavior  is  isotropic,  the  elastic  moduli  can  be  expressed  as 

CJh  =  +  w  +  (k2g )  5ijdkl  (17) 

where  G  and  I<  are  the  elastic  shear  and  bulk  moduli,  respectively.  Therefore,  to  update  stresses,  the  plastic  strain  increments 
need  to  be  determined.  The  following  outlines  the  procedure  for  computing  A ejj. 

The  yield  condition  and  the  flow  rule  are  written  as 

F(l\+At,J2+AtJ3At)-S(ept+At)=  0  (18) 

and 


Eqs.  (18)  and  (19),  when  considered  together  with  Eq.  (15),  result  in  10  equations  for  A2  and  nine  components  of  A ejj, 
among  which  seven  equations  are  independent  due  to  the  symmetry  of  Aejj.  If  the  state  variable  W  is  updated  and  thus 
d(ef+At)  is  known,  these  equations  can  be  solved  iteratively  for  A .e?.  and  A2  using  the  Newton-Raphson  method.  The  iterative 
process  follows  these  steps:  (1)  assume  Aejj  =  0  and  use  Eq.  (15)  to  estimate  a\PAt;  (2)  use  Eqs.  (18)  and  (19)  to  solve  for  A ejj; 
(3)  update  stresses  using  Eq.  (15);  (4)  repeat  steps  (2)-(3)  until  convergence  conditions  are  satisfied. 

To  update  sp ,  consider  the  hardening  law  and  the  evolution  equation  for  sp 


0-  =  S(spt+At)  =  a(ept  +  Asp) 

(20) 

gAIp  =  G\PAtA8?j 

(21) 

At  each  iteration  of  A ejj,  A sp  and  <j(spt+At)  can  be  obtained  by  solving  Eqs.  (20)  and  (21 )  iteratively  using  the  Newton-Raph¬ 
son  method. 


2.2.2.  The  consistent  tangent  moduli 

Simo  and  Taylor  (1985)  showed  that  use  of  the  consistent  tangent  moduli  significantly  improves  the  convergence  char¬ 
acteristics  of  the  overall  equilibrium  iterations.  The  consistent  tangent  stiffness  corresponding  to  the  backward  Euler  inte¬ 
gration  can  be  obtained  by  linearization  of  Eq.  (15). 


Jiikl  (dekl 


(22) 


Since  all  quantities  in  calculating  Jp*  are  referred  to  time  t  +  At,  the  superscript  t  +  At  will  be  dropped  from  hereafter. 
Eq.  (15)  can  be  rewritten  as 
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°ij  =  Cijki4i  =  Gjki(Ski  -  Spkl)  =  C|M{gw  -  (eg, )‘  -  AgP  } 

Differentiating  (23)  results  in 

dOij  =  CjjudEki  -  C|wa(AfiP ) 

The  relationship  between  day  and  d(A£^d)  can  be  found  as  follows. 

Let  h  be  the  hardening  modulus,  i.e.,  h  =  Aa/AeP,  then  Eq.  (21)  can  be  rewritten  as 


Act  = 


hVijA  eg 


By  differentiating  (25)  and  simplifying  the  resulted  relation,  the  following  equation  can  be  obtained 

ha  a  a 


da  =  — 


haAs^ 


-da«  +  — 


-d 


H) 


a2  +  hamnA£pmn  lJ  a2  +  hamriA£pmn 

Differentiating  Eq.  (18)  and  combining  the  result  with  Eq.  (26)  give 

haau 


a2  +  hi JmnAspmn 


-d 


K) 


mh.J2J3)  haAtf 

a2  +  hamnAspmr  1  l] 


da , 


Eliminating  AA  from  the  nine  equations  given  by  (19)  results  in  eight  equations,  such  as 


(23) 

(24) 

(25) 

(26) 

(27) 


(28) 


Differentiating  (28)  leads  to 

-  (S',<Ae«  - 

ds><As!'>  -  {&)d^  - 


A  P  d2G 

822  dOndOij 

a  p  d2° 

A  22  dOjrfOij 


-M, 


d2C 

d<J22d(Jij 

A  P  926 

^802280, 


dan 


dan 


(29) 


Eqs.  (27)  and  (29)  provide  nine  equations  between  and  day,  which  can  be  summarized  as 

IC9(A£p)  =  Dcto  (30) 

where  Asp  =  {Ae^ ,  Ae^1 ,  Ae^ ,  •  •  • ,  Agf3 }T,  a  =  {au ,  a2\ ,  cr3i ,  •  •  • ,  a33}T,  K  is  the  coefficient  matrix  of  d(Asp)  and  D  is  the  coef¬ 
ficient  matrix  of  da. 

From  Eqs.  (24)  and  (30)  we  can  obtain 

d{  Asp)  =  (K  +  DCy'DCide)  (31) 

where  Ce  is  a  9  x  9  matrix  representing  the  elasticity  tensor  C|w.  Finally,  the  consistent  tangent  matrix  can  be  obtained  by 
substituting  (31)  into  (24) 

J  =  9<T/a«  =  Ce-Ce(K  +  DCe)'1DC£'  (32) 


2.3.  A  modified  Gurson-Tvergaard-Needleman  model 

Ductile  fracture  of  many  structural  materials  is  a  result  of  void  nucleation,  growth  and  coalescence.  The  constitutive 
description  of  this  mechanism  has  received  a  great  deal  of  attention  in  the  past  30  years,  which  leads  to  various  forms  of 
porous  material  models  being  developed  to  describe  void  growth  and  the  associated  macroscopic  softening.  One  of  the  most 
famous  porous  plasticity  models  was  due  to  Gurson  (1977)  with  modifications  by  Tvergaard  and  Needleman  (Tvergaard, 
1981,  1982;  Tvergaard  and  Needleman,  1984).  In  the  Gurson-Tvergaard-Needleman  model,  an  extra  internal  variable,  the 
void  volume  fraction  (/),  is  introduced  to  capture  the  growth  of  cavities  and  its  effect  on  material  behavior.  It  is  important 
to  notice  that  the  Gurson-Tvergaard-Needleman  model  reduces  to  that  for  J2-flow  theory  of  plasticity  with  isotropic  hard¬ 
ening  in  the  absence  of  voids  (f=  0). 
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Within  the  framework  of  the  l\-J2-J3  plasticity  theory  outlined  in  Section  2.1,  the  yield  function  and  flow  potential  of  the 
Gurson-Tvergaard-Needleman  model  can  be  modified  as 

*  =  (£) 2  +  2<3i/cosh  -  1  -  q2J2  =  0  (33) 

and 


W  = 


+  2q1/cosh 


q?/2  =  o 


(34) 


where  F and  G  take  the  forms  of  Eqs.  (6)  and  (9),  i.e.,  F  =  C\  (a^  +  27/|  +  bif^j  1  and  G  =  c2(a2I \  +  27/2  +  bj^j  '  ,  /  J2  and 

J3  are  invariants  of  the  macroscopic  stress,  /  is  the  current  void  volume  fraction,  a  is  the  current  yield  stress  of  the  matrix 
material,  and  qx  and  q2  are  parameters  introduced  by  Tvergaard  (1981,  1982)  to  account  for  void  interaction  and  matrix 
strain  hardening.  If  a\  =  b\  =  a2  =  b2  =  0,  Eqs.  (33)  and  (34)  degenerate  to  the  original  Gurson-Tvergaard-Needleman  model. 

Porous  material  models  contain  an  additional  state  variable,/.  For  the  modified  Gurson-Tvergaard-Needleman  model,  the 
evolution  equation  for /can  be  obtained  by  considering  the  rate  of  the  net  volume  change 


/=(!-/) 


4-3^2/,^ 


(35) 


By  enforcing  equality  between  the  rates  of  macroscopic  plastic  work  and  the  matrix  plastic  dissipation,  the  matrix  yield 
stress,  (T,  and  the  matrix  plastic  strain  rate,  P,  are  coupled  through 


(1  -f)aW  = 


(36) 


where  the  matrix  material  follows  a  prescribed  hardening  function  a(sp). 

Numerical  implementation  of  the  modified  Gurson-Tvergaard-Needleman  model  follows  a  similar  procedure  as  outlined 
in  Section  2.2. 


3.  Numerical  examples 

3.1.  Modeling  the  plastic  behavior  of  an  aluminum  5083  alloy 

In  this  section  we  apply  the  t\-J2-]3  plasticity  model  to  study  the  plastic  response  of  an  aluminum  alloy  5083-H116, 
which  was  cold  worked  to  achieve  its  strength.  All  specimens  are  machined  from  a  25  mm  thick  plate,  with  tensile  axes  ori¬ 
ented  transversely  to  the  rolling  direction.  Round  specimens  were  turned  with  low  stress  machining  procedures  and  rectan¬ 
gular  specimens  were  electro-discharge  machined.  All  tests  are  performed  at  room  temperature  and  are  considered  to  be 
quasi-static.  The  test  matrix  includes  smooth  and  notched  round  tensile  bars,  grooved  plane  strain  specimens  and  the  Lind- 
holm-type  specimen  (Lindholm  et  al.,  1980)  subjected  to  different  tension-torsion  ratios.  Fig.  1  shows  a  smooth  round  bar,  a 
notched  round  bar,  a  grooved  plane  strain  specimen  and  a  torsion  specimen. 

Two  non-dimensional  parameters  defined  as  T  =  f  /  (3V3f2/2^j  and  £  =  3V3J3/ together  with  the  equivalent  stress 

<je,  are  often  used  to  characterize  the  stress  state,  where  T  is  the  stress  triaxiality  ratio  and  £  can  be  related  to  the  Lode  param¬ 
eter  (Wierzbicki  and  Xue,  2005;  Xue,  2007;  Kim  et  al.,  2004,  2007;  Gao  et  al.,  2005,  2009,  2010).  For  the  smooth  and  notched 
round  tensile  bars,  the  Rvalue  at  the  center  of  the  specimen  is  one  while  the  T-value  varies  with  the  notch  radius.  The  smal¬ 
ler  the  notch  radius  is,  the  higher  the  stress  triaxiality.  On  the  other  hand,  the  Rvalue  at  the  center  of  the  grooved  plane 

Fig.  1.  Sketches  of  a  smooth  round  bar,  a  notched  round  bar,  a  grooved  plane  strain  specimen  and  a  torsion  specimen. 
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strain  specimens  is  zero  and  the  T-value  varies  with  the  groove  radius.  The  Lindholm-type  specimen  can  achieve  various 
combinations  of  T  and  £  by  varying  the  ratio  between  the  applied  tensile  displacement  and  the  applied  twist  angle. 

The  diameter  of  the  gage  section  of  the  smooth  round  bar  is  6.35  mm.  All  the  notched  round  bars  have  the  same  diameters 
of  15.2  mm  in  the  smooth  sections  and  7.6  mm  at  their  notched  cross  section.  Three  notch  radii,  1.27,  2.54  and  6.35  mm,  are 
considered  for  specimens  D,  B  and  E,  respectively  (Table  1 ).  The  overall  dimensions  for  the  circular-grooved  plane  strain 
specimens  are  203.2  x  31.8  x  6.1  mm  and  the  thinnest  cross  section  is  2  mm.  Three  different  groove  radii,  2.03,  5.08  and 
16.26  mm,  are  considered  for  specimens  F,  G  and  H,  respectively  (Table  2).  The  torsion  specimens  are  hollow  cylinders  hav¬ 
ing  an  inner  diameter  of  13.1  mm  and  outer  diameter  of  23.8  mm.  The  gage  section  has  a  length  of  2.54  mm  and  wall  thick¬ 
ness  of  0.75  mm.  The  torsion-tension  specimens  have  a  slightly  larger  outer  diameter  at  both  ends  (25.4  mm)  and  other 
dimensions  are  the  same  as  the  torsion  specimen.  Table  3  lists  the  ratio  between  the  applied  tensile  displacement  and  ap¬ 
plied  twist  angle  for  each  of  the  tension-torsion  specimens.  Details  of  specimen  drawings  are  given  in  Gao  et  al.  (2009). 

The  specimens  in  the  test  matrix  of  this  study  generate  a  wide  range  of  stress  states.  Moreover,  the  stress-state  of  a  mate¬ 
rial  point  in  the  test  specimens  evolves  as  plastic  deformation  increases.  Figs.  2-5  show  the  variation  of  T  =  A  /  (3  V3j]/2^j  and 

£  =  3\/3J3/ (2jl/2^)  with  loading  history  (measured  by  the  equivalent  plastic  strain,  sp)  in  the  element  at  the  specimen  center 

for  the  smooth  round  bar,  the  E-notch  round  bar,  the  G-groove  plane  strain  specimen  and  the  TT-16  torsion-tension  spec¬ 
imen.  For  the  smooth  round  bar,  £  remains  at  1  during  the  entire  loading  history  while  T  increases  from  1/3  to  about  0.45 
before  specimen  fractures.  For  the  E-notch  specimen,  £  remains  at  about  1  and  T  increases  from  0.71  to  0.8.  For  the  G-groove 
specimen,  £  quickly  decreases  to  0  and  remains  at  this  level  as  plastic  deformation  increases  while  T  increases  from  0.51  to 
0.8  before  failure  occurs.  For  the  tension-torsion  specimen  TT-16,  £  increases  from  0.4  to  0.94  and  T  increases  from  0.1  to 
0.26.  For  the  notched  and  grooved  specimens,  changing  notch  (groove)  radius  changes  the  level  of  T  in  the  specimen.  Con¬ 
sidering  the  entire  loading  history  of  each  specimen  and  the  three  different  notch  (groove)  radii,  the  range  of  T  experienced 
by  the  center  element  is  0.71  <  T<  1.6  for  the  notched  round  bar  tests  and  0.46  <  T  ^  0.97  for  the  plane  strain  tests.  For  the 
tension-torsion  tests,  varying  the  ratio  between  the  applied  tensile  displacement  and  the  applied  twist  angle  results  in  dif¬ 
ferent  combinations  of  the  £  and  T  histories.  Considering  the  entire  loading  history  and  the  different  tensile  displacement/ 


Table  1 

Notch  radii  of  the  notched  round  bars. 


Specimen 

B-notch 

D-notch 

E-notch 

Notch  radius  (mm) 

2.54 

1.27 

6.35 

Table  2 

Groove  radii  of  the  plane  strain  specimens. 


Specimen 

F-groove 

G-groove 

H-groove 

Groove  radius  (mm) 

2.03 

5.08 

16.26 

Table  3 

Ratios  of  the  applied  tensile  displacement  and  applied  twist  angle  used  in  the  tension-torsion  tests. 


Specimen 

TT-17 

TT-19 

TT-13 

TT-16 

TT-15 

TT-14 

tensile  displacement/twist  angle  (mm/radian) 

0.23 

0.69 

0.92 

1.16 

1.38 

2.54 
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Fig.  2.  Variation  of  T  =  G/  (3V3j]/2^  and  £  =  3 v/3/s /  with  plastic  deformation  in  the  center  element  of  the  smooth  round  bar. 
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Fig.  3.  Variation  of  T  =  G/ (3V3Jl/2^  and  £  =  3  v/3/3 /  (^2 /|/2 ^  with  plastic  deformation  in  the  center  element  of  the  E-notch  specimen. 


Fig.  4.  Variation  of  T  =  I\/(3 V3jl/2^)  and  £  =  3V3J3/ (2jl/2>)  with  plastic  deformation  in  the  center  element  of  the  G-groove  specimen. 


Fig.  5.  Variation  of  T  =  G  /(3y/3Jl/2^  and  £  =  3 V3J3 /(2Jl/2>)  with  plastic  deformation  in  the  center  element  of  the  tension-torsion  specimen  TT-16. 


twist  angle  ratios  listed  in  Table  3,  the  range  of  T  and  £  experienced  by  the  tension-torsion  specimens  in  the  test  matrix  are 
in  the  range  of  0.018  <  0.365  and  0.08  ^  £  <  0.99. 

Numerical  analyses  of  the  specimens  are  carried  out  using  the  general  purpose  finite  element  software  ABAQUS  (SIMULIA, 
2008),  which  employs  an  updated  Lagrangian  formulation  to  handle  finite  deformation.  The  material  models  outlined  in  pre¬ 
vious  sections  are  implemented  in  ABAQUS  via  user  defined  subroutines.  For  round  tensile  bars,  axisymmetric  conditions  are 
considered  and  the  4-node,  axisymmetric  solid  elements  with  reduced  integration  (CAX4R)  are  used.  For  torsion-tension 
specimens,  an  additional  degree  of  freedom  needs  to  be  added  to  the  axisymmetric  element  to  handle  the  twist  and  the 
CGAX4R  element  in  ABAQUS  is  developed  for  this  purpose.  For  grooved  plane  strain  specimens,  the  3D  8-node  brick  elements 
with  reduced  integration  (C3D8R)  are  used.  Usually  the  symmetry  conditions  allows  for  only  1/4  or  1/8  of  the  specimen 
being  modeled.  Fig.  6  shows  two  typical  finite  element  meshes.  A  typical  axisymmetric  model  has  700  elements  and  a  typical 
1  /8-symmetric  3D  model  has  20,000  elements.  Since  reduced  integration  is  used,  the  hourglass  control  is  employed  by  using 
HOURGLASS  STIFFNESS  in  the  ABAQUS  input  files  to  provide  increased  resistance  to  hourglassing. 

Fig.  7  shows  the  true  (von  Mises)  stress  vs.  plastic  strain  curves  obtained  using  the  smooth  tensile  bar  data  and  the  pure 
torsion  test  data.  Gao  et  al.  (2009)  provides  the  details  of  how  these  two  curves  were  obtained  from  the  measured  load-dis¬ 
placement  and  torque-twist  angle  curves,  respectively.  The  two  curves  clearly  show  significant  difference,  suggesting  that 
the  stress-state  has  a  strong  effect  on  the  plastic  response  of  this  material. 
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Fig.  6.  Typical  finite  element  meshes:  (a)  an  axisymmetric  model  for  a  notched  round  bar,  (b)  a  1  /8-symmetric  model  for  a  grooved  plane  strain  specimen. 


Plastic  Strain 


Fig.  7.  Comparison  of  the  true  (von  Mises)  stress  vs.  true  plastic  strain  curves  obtained  using  the  smooth  tensile  bar  data  and  the  torsion  test  data. 

Now  consider  the  l\-]2~h  plasticity  model  presented  in  previous  sections  and  conduct  finite  element  analyses  of  these 
two  specimens.  The  uniaxial  tension  stress-strain  curve  can  be  described  by  a  power-law  hardening  relation 

£  =  f  when  g  <  a0 

/  \VN  (37) 

£  =  T  (%)  when  °  >  ^0 

with  E  =  68.4  GPa,  g0  =  198.6  MPa,  v  =  0.3,  and  N=  0.155.  The  plasticity  model  requires  a  relation  of  the  current  yield  stress  as 
a  function  of  the  effective  plastic  strain.  This  is  achieved  by  using  a  UHARD  subroutine.  Using  the  uniaxial  tension  stress- 
strain  curve  and  adjusting  material  parameters  (cq.ibi)  and  (a2,  b2),  it  is  found  that  a  set  of  parameters,  ai  =  a2  =  0, 
b\  =  -60.75  and  b2  =  -25,  result  in  a  best  match  between  the  numerically  predicted  and  experimentally  measured  torque 
vs.  twist  angle  response  of  the  pure  torsion  test.  This  indicates  no  Ji  effect  but  significant  J3  effect  on  the  plastic  response. 
Fig.  8  compares  the  numerical  predictions  and  the  experimental  measurements  for  the  smooth  tensile  specimen  and  the 
pure  torsion  specimen,  respectively.  The  tensile  tests  have  seven  specimens  and  the  torsion  tests  have  10  specimens.  The 
experimental  data  are  represented  by  thinner  lines  while  the  finite  element  result  is  represented  by  a  thicker  line.  In  both 
cases,  excellent  comparison  between  the  numerical  prediction  and  the  experimental  measurements  is  observed. 

With  the  parameters  of  the  plasticity  model  being  calibrated,  the  next  question  needs  to  be  answered  is  whether  this 
model  correctly  predicts  the  material  response  under  complex  stress  states.  To  this  end,  notched  round  bars  with  different 
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Fig.  8.  Comparison  of  the  numerical  and  experimental  load  vs.  displacement  curves  for  (a)  the  smooth  tensile  specimen  and  (b)  the  pure  torsion  specimen. 


notch  radii,  plane  strain  specimens  with  different  groove  radii  and  torsion  specimens  subjected  to  different  tension-torsion 
ratios  are  tested  and  analyzed  and  the  numerical  predictions  are  compared  with  the  experimental  records.  For  notched 
round  bars  and  grooved  plane  strain  specimens,  the  applied  tensile  force  vs.  extensometer  gage  displacement  response  is 
monitored.  For  tension-torsion  specimens,  both  applied  axial  force  vs.  axial  displacement  and  applied  torque  vs.  twist  angle 
responses  are  monitored.  Fig.  9  shows  typical  comparisons  -  again,  model  predictions  agree  with  experimental  measure¬ 
ments  very  well.  Similar  comparisons  between  experimental  results  and  numerical  predictions  are  observed  for  all  other 
specimens  in  the  test  matrix. 

To  further  justify  the  proposed  plasticity  model,  Fig.  10a  and  b  compares  the  numerical  predictions  using  the  classical  J2- 
flow  theory  and  the  proposed  /1-J2-J3  model  for  the  tension-torsion  test  TT-16  with  experimental  results.  In  these  figures, 
we  also  included  the  comparison  between  the  non-associated  flow  rule  (Non-AFR)  and  the  associated  flow  rule  (AFR).  As  can 
be  seen  from  the  plots,  only  the  proposed  l\-J2-h  model  with  the  calibrated  parameters  (non-associated  flow  rule)  leads  to 
good  agreement  with  experimental  records  for  both  torque  vs.  twist  angle  and  axial  force  vs.  axial  displacement.  The  asso- 
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Fig.  9.  Comparisons  of  the  numerical  predictions  and  experimental  records:  (a)  notched  round  bar  (E-Notch);  (b)  plane  strain  specimen  (G-Groove);  (c) 
torque  vs.  twist  angle  response  for  torsion-tension  test  (TT-16);  (d)  axial  force  vs.  axial  displacement  response  for  torsion-tension  test  (TT-16). 
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Fig.  10.  Comparisons  of  the  numerical  predictions  using  the  classical  J2-flow  theory  and  the  proposed  U-J2-J3  model  for  the  tension-torsion  test  TT-16  with 
experimental  results:  (a)  torque  vs.  twist  angle;  (b)  axial  force  vs.  axial  displacement. 


ciated  flow  rule  gives  acceptable  torque  vs.  twist  angle  prediction  but  unsatisfying  axial  force  vs.  axial  displacement 
prediction. 

Fig.  1 1  shows  the  difference  between  the  locus  of  the  plastic  flow  potential  (G)  and  the  locus  of  the  yield  function  (F)  given 
by  the  /1-J2-J3  model  with  aA  =  a2  =  0,  bA  =  -60.75  and  b2  =  -25  on  the  G\-o2  plane,  where  the  locus  of  the  flow  potential  is 
represented  by  the  solid  line  and  the  locus  of  yield  function  is  represented  by  the  dashed  line. 

In  summary,  the  plasticity  behavior  of  the  5083  aluminum  alloy  considered  in  this  study  depends  on  the  stress  state  and 
can  be  described  by  an  U-J2-J3  plasticity  model  described  in  Section  2.  The  calibrated  model  constants  are  a.\  =  a2  =  0, 
hi  =  -60.75  and  b2  =  -25,  indicating  no  U  effect  but  significant  J3  effect  on  the  plastic  response  of  this  material. 


3.2.  Numerical  examples  using  the  Gurson-Tvergaard-Needleman  model 

In  this  section  we  conduct  a  series  of  parametric  studies  to  illustrate  the  effect  of  the  modified  Gurson-Tvergaard-Nee¬ 
dleman  model  on  predicted  material  response.  Fig.  12  shows  a  cubic  element  with  dimensions  of  D0  x  D0  x  D0.  Displacement 
boundary  conditions  are  imposed  on  the  element  surfaces  such  that  the  macroscopic  stress  ratios,  pi  =  a11/a22  and  p2  =  cr33/ 
(722,  are  kept  constants  during  the  entire  deformation  history  (which  results  in  a  constant  stress  triaxiality  ratio  and  Lode 
angle).  Faleskog  et  al.  (1998)  and  Kim  et  al.  (2004)  provide  details  of  how  to  prescribe  this  kind  of  boundary  conditions. 
The  material  parameters  used  in  the  analyses  are  E  =  68.4  GPa,  0  =  0.3,  <r0  =  207  MPa,  N  =  0.125  and  f0  =  0.003,  where  F,  0, 
(70,  n  and/o  represent  the  Young’s  modulus,  the  Poisson’s  ratio,  the  yield  stress,  the  strain  hardening  exponent  and  the  initial 
void  volume  fraction,  respectively.  The  qA  and  q2  parameters  in  the  Gurson-Tvergaard-Needleman  model  are  taken  as 
q\  —  1.5  and  q2  =  \. 


Fig.  11.  Comparison  between  the  locus  of  the  plastic  flow  potential  (G)  and  the  locus  of  the  yield  function  (F)  given  by  the  U-J2-J3  model  with  a1  =  a2  =  0, 
hi  =  -60.75  and  b2  =  -25  on  the  0\-o2  plane. 
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For  the  first  set  of  analyses,  the  boundary  conditions  are  imposed  such  that  p ^  =  0.268  and  p2  =  0.634,  corresponding  to  a 
stress  triaxiality  ratio  of  1  and  Lode  angle  of  0°.  When  a^=a2  =  b^  =  b2  =  0,  the  analysis  results  using  our  user  subroutine  are 
the  same  as  those  obtained  using  the  original  Gurson-Tvergaard-Needleman  model  implemented  in  ABAQUS,  serving  as  a 
check  of  our  numerical  implementation. 

Fig.  13  illustrates  the  effect  of  A,  where  the  dotted  lines  represent  the  results  obtained  using  aA  =  a2  =  b\  =  b2  =  0  and  the 
solid  lines  represent  the  results  obtained  using  ci\  =  a2  =  6  x  10-4  and  b\  =  b2  =  0.  At  the  same  applied  displacement  in  the  x2- 
direction  (u2)»  both  o22  and  the  void  growth  rate  are  lower  when  the  effect  of  A  is  taken  into  account. 


a  b 

^22  /  ^0  f 


Fig.  13.  Comparison  of  the  <j22I<7q  vs.  u2/D0  response  and /vs.  u2ID0  response  predicted  using  ai  =  a2  =  hi  =  b2  =  0  (dotted  lines)  and  a^  =  a2  =  6  x  10  4  and 
b^=b2  =  0  (solid  lines). 


a  b 
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Fig.  14.  Comparison  of  the  a22la0  vs.  u2/D0  response  and  /  vs.  u2/D0  response  predicted  using  a-i  =  a2  =  b-i  =  b2  =  0  (dotted  lines)  and  ai  =  a2  =  0  and 
b-l=b2  =  -60.75  (solid  lines). 
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a  b 
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Fig.  15.  Comparison  of  the  <722/00  vs.  u2/D0  response  and /vs.  u2/D0  response  predicted  using  ai  =  a2  =  6xl0  4  and  b^  =  b2  =  0  (associated  flow  rule;  dotted 
lines)  and  at=6x  10-4  and  a2  =  b^  =  b2  =  0  (non-associated  flow  rule;  solid  lines). 


Fig.  14  demonstrates  the  effect  of  J3  by  setting  non-zero  values  for  bA  and  b2.  Here  the  dotted  lines  show  the  results  of 
a.\  =  a2  =  b\  =  b2  =  0  and  the  solid  lines  are  obtained  by  using  aA  =  a2  =  0  and  b^  =  b2  =  -60.75.  The  analysis  results  show  that 
negative  values  of  b\  and  b2  lead  to  lower  value  of  cr22  and  slower  void  growth  rate. 

The  results  shown  in  Figs.  13  and  14  are  obtained  using  the  associated  flow  rule.  Figs.  15  and  16  compare  the  results  of  the 
non-associated  flow  rule  with  those  of  the  associated  flow  rule.  As  is  shown  in  Fig.  15,  where  the  dotted  lines  are  the  results 
of  an  associated  flow  rule  (ai  =  a2  =  6  x  10-4  and  b\  =  b2  =  0)  and  the  solid  lines  are  the  results  of  a  non-associated  flow  rule 
(ai  =  6  x  10-4  and  a2  =  b\  =  b2  =  0),  increasing  a2  from  0  to  6  x  10-4  leads  to  an  increase  of  o22  and  decrease  of  void  growth 
rate. 

Fig.  16  compares  the  results  when  b2  takes  a  different  value  than  b\.  In  these  analyses,  aA  and  a2  are  taken  as  zero  and  the 
results  using  b^  =  b2  =  -60.75  (associated  flow  rule)  are  denoted  by  the  dotted  lines  while  the  results  using  bA  =  -60.75  and 
b2  =  0  (non-associated  flow  rule)  are  denoted  by  the  solid  lines.  Varying  b2  from  0  to  -60.75  results  in  negligible  effect  on  o22 
but  a  higher  void  growth  rate. 

In  the  analyses  presented  above,  the  boundary  conditions  are  imposed  to  keep  pi  =  0.268  and  p2  =  0.634.  The  stress-state 
experienced  by  the  material  can  be  altered  by  changing  the  boundary  conditions  so  that  different  values  of  pi  and  p2  are 
achieved.  For  example,  p  1  =  0.4  and  p2  =  0.4  corresponds  to  a  stress  triaxiality  ratio  of  1  and  Lode  angle  of  -30°.  Fig.  17  illus¬ 
trates  how  this  change  of  stress  state  (Lode  angle  changes  from  0°  to  -30°  while  stress  triaxiality  ratio  remains  1 )  affects  o22 
and  the  void  growth  rate,  where  a\  =  a2  =  6  x  10-4  and  b^  =  b2  =  -60.75  are  used  in  the  calculations.  In  Fig.  17,  the  dotted 
lines  represent  the  results  of  pi  =  0.268  and  p2  =  0.634  (Lode  angle  equal  to  0°)  and  the  solid  lines  represent  the  results  of 
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^22  /  ^0 
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Fig.  16.  Comparison  of  the  022/00  vs.  u2/D0  response  and/vs.  u2/D0  response  predicted  using  bA  =  b2  =  -60.75  (dotted  lines)  and  b\  =  -60.75  and  b2  =  0  (solid 
lines),  where  and  a2  are  taken  as  zero. 
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Fig.  17.  Comparison  of  the  <722/00  vs.  u2/D0  response  and /vs.  u2ID0  response  when  pi  =  0.268  and  p2  =  0.634  (dotted  lines)  with  those  when  pi  =  0.4  and 
p2  =  0.4  (solid  lines),  where  a1=a2  =  6x  10-4  and  b\=b2  =  -60.75. 


Pi  =  0.4  and  p2  =  0.4  (Lode  angle  equal  to  -30°).  The  void  growth  rate  is  slightly  higher  and  cr22  becomes  noticeably  larger 
when  the  Lode  angle  changes  from  0°  to  -30°. 

In  summary,  the  dependence  of  the  matrix  plasticity  behavior  on  U  and  J3  results  in  significant  changes  in  void  growth 
and  the  macroscopic  stress  vs.  deformation  response  of  porous  materials.  The  modified  Gurson-Tvergaard-Needleman  mod¬ 
el  provides  a  means  to  describe  these  changes. 


4.  Concluding  remarks 

In  this  paper,  we  describe  a  plasticity  model  for  isotropic  materials,  which  is  dependent  on  the  second  and  third  invariants 
of  the  stress  deviator  as  well  as  the  hydrostatic  stress,  and  present  its  finite  element  implementation  including  integration  of 
the  constitutive  equations  using  the  backward  Euler  method  and  the  formulation  of  the  consistent  tangent  moduli.  As  an 
application,  this  model  is  calibrated  and  verified  for  a  5083  aluminum  alloy.  The  model  captures  the  strong  stress-state  effect 
on  the  plastic  behavior  of  this  material  and  accurately  predicts  the  plastic  responses  of  specimens  experiencing  a  wide  range 
of  stress  states.  Furthermore,  the  Gurson-Tvergaard-Needleman  porous  plasticity  model,  which  is  widely  used  to  simulate 
the  void  growth  process  of  ductile  fracture,  is  extended  to  include  the  effects  of  the  third  invariant  of  the  stress  deviator  and 
the  hydrostatic  stress  on  the  matrix  material.  A  series  of  parametric  studies  illustrate  the  effects  of  model  parameters  on  the 
predicted  material  response.  This  simple  modification  enriches  the  Gurson-Tvergaard-Needleman  model  and  expands  its 
applicability  as  a  micromechanical  model  for  ductile  fracture. 
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